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ABSTRACT 

Observations in the past decade have revealed extrasolar planets with a wide range of orbital semi- 
major axes and eccentricities. Based on the present understanding of planet formation via core accre- 
tion and oligarchic growth, we expect that giant planets often form in closely packed configurations. 
While the protoplanets are embedded in a protoplanetary gas disk, dissipation can prevent eccen- 
tricity growth and suppress instabilities from becoming manifest. However, once the disk dissipates, 
eccentricities can grow rapidly, leading to close encounters between planets. Strong planet-planet 
gravitational scattering could produce both high eccentricities and, after tidal circularization, very 
short-period planets, as observed in the exoplanet population. We present new results for this sce- 
nario based on extensive dynamical integrations of systems containing three giant planets, both with 
and without residual gas disks. We assign the initial planetary masses and orbits in a realistic manner 
following the core accretion model of planet formation. We show that, with realistic initial conditions, 
planet -planet scattering can reproduce quite well the observed eccentricity distribution. Our results 
also make testable predictions for the orbital inclinations of short-period giant planets formed via 
strong planet scattering followed by tidal circularization. 

Subject headings: methods: n-body simulations, methods: numerical, (stars:) planetary systems, 
(stars:) planetary systems: protoplanetary disks, planetary systems: formation — 
celestial mechanics 



1. INTRODUCTION 

The study of extrasolar planets and their properties 
has become a very exciting area of research over the 
past decade. Since the detec tion of the planet 5 1 Peg b, 
more than 200 new p lanets (jButler et al.l 120(361 see also 
http://exoplanet.eu/) have been detected and the large 
sky surveys planned for the near future can potentially 
detect many more. These detections have raised many 
questions about the formation and dynamical evolution 
of planetary systems. The extrasolar planet population 
covers a much greater portion of the semimajor axis and 
eccentricity plane than was expected based on the plan- 
ets in our solar system ([Lissa ucr 1995, Fig. [I}. The pres- 
ence of many giant planets in highly eccentric orbits or 
in very short-period orbits (the "hot Jupiters" ) is partic- 
ularly puzzling. 

Different scenarios have been proposed to explain the 
high eccentricities. The presence of a distant compan- 
ion in a highly inclined orbit can increase the eccentric- 
ities o f the planets around a star through Ko zai oscilla- 
tions (|Mazeh et al.lll997tlHolman et al.lll997ft . However, 
this alone cannot explain the ob served eccentricity dis- 
tribution ( Ta keda fe Rasiol [2005) . Interaction with the 
protoplanetary gas disk could either excite or damp the 
eccentricities depending on the properties of the disk 
and the orbits. However, the com bined effects typi- 
cally result in eccentric i ty damping (|Artvmowic3 1 1992 1: 
Papaloizou fc TerquerrJ 120011: IGoldreich fc Saril l20ul 
Ogilvie fc Lubo w 2003). Migration of two planets and 



trapping in a mean-motion resonance (MMR) can also 
pump up the eccentricities efficiently, but this mech- 
anism requires strong damping at the end or termi- 
nation of migratio n right after trapping in resonance 
( Lee fc Peald l2002p or el se it leads to planet sca t terin 
jSandor fc Kiev! l2006ft . IZakamska fc Tremaind (|200 
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proposed inward propagation of eccentricity after the 
outer planets are excited to high eccentricities following a 
close encounter with a passing star. Using typical values 
for such interactions with field stars in the solar neighbor- 
hood, however, they do not ge t very high eccentricities. 
iPapaloizou fc Terqueml (|200lD : iTerquem fc Papaloizou! 
(12002ft : lBlackl (j!997ft propose a very different formation 
scenario for planets from protostellar collapse in which 
both hot Jupiters and eccentric planets at higher semi- 
major axes are formed naturally; this scenario, however, 
cannot form sub- Jupiter-mass planets. 

In this paper, we explore another promising way 
to create high eccentricities: strong gravitational scat- 
tering between planets in a multi -planet system un- 
dergoing dynamical instability dRasio fc Fordl 119961 : 
IWeidenschilling fc Marzarll 119961 : ILin fc Idal 11997ft . Ac- 
cording to the model of oligarchic growth, planetesi- 
mals form in a nearly maximally packed configuration 
in the protoplanetary disk, followed by gas accretion 
(IGoldreich et al.ll20ui Hda fc Linl [2004bT iKokubo fc Idal 
12002ft . Once the disk dissipates, mutual planetary per- 
turbations ( "viscous stirring" ) of the planetesimals will 
lead to eccentricity growth, orbit crossing, and even- 
tually ^_da3e_encouirte^ in th e 
disk (|Ford fc Chiand [20071 : Eevison fc Morbidellil 12007ft . 
While planetary systems with more than two planets 
can not be provably stable, they can remain stable for 
very l ong timescales depending on their initial separa- 
tions (Chamb ers et al.lll996l : iMarzari fc Weidensc hilliiig 
2002). A sufficiently massive disk can prevent interacting 
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Fig. 1. — Semi-major axis vs eccentricity for the detected planets. 
The sizes of the points are proportional to the minimum masses 
(msini) of the corresponding planets. The size of a Jupiter mass 
planet is shown at the left top corner for reference. The stars 
represent the four giant planets in our solar system (for these the 
sizes do not indicate their mass). The open circles show the planets 
detected by radial velocity surveys. The triangles show planets 
detected by micro-lensing or direct imaging. Planets with poorly 
constrained eccentricities are plotted above e = 1. A horizontal 
line is drawn at e = 1 to guide the eye. Note the logarithmic scale 
for the semi-major axis. 

planets from acquiring large eccentricities and developing 
crossing orbits. However, once the gas disk is sufficiently 
dissipated, and the planetesimal disk depleted, the eccen- 
tricities of the planets can grow to high values, possibly 
leading to strong planet-planet scattering and a phase 
of chaotic evolution tha t dramatically alters the orbital 
structure of the system (|Lin fc Idalll997tlFord fe ChianeJ 
120071 : iLevison fc Morbideliill2007h . 

The detection of close-in planets with orbital periods 
as short as ~ Id, the so-called "hot Jupiters" (and, 
more recently, hot Neptunes and super-Earths) , was an- 
other major surprise. Giant planets are most likely to 
form at much larger separations, beyond the ice line of 
the star where there can be enhanced dust production 
(jKokubo fc Idall200llida fc Linll2004br ). It is widely be- 
lieved that the giant planets form beyond the ice line and 
then migrate inwards to form the hot Jupiters we observe 
today. Different stopping mechanisms of inward migra- 
tion have been proposed to explain the hot Jupiters, 
but it is unclear why they pile up at just a few solar 
radii around the star, rather than continue migrating and 
eventually accrete onto the star. 

Strong gravitational scattering between planets in a 
multi-planet system may provide another way to cre- 
ate these close-in planets ?Rasio fc Fordl fT996l ) . A few 
of the planets scattered into very highly eccentric orbits 
could have sufficiently small periastron distances that 
tidal circularization takes place, giving rise to the hot 
Jupiters. The currently observed edge in the mass-period 
diagram is very nearly at the ideal circularization radius 
twice the Roche lim i t) , providing suppo rt for this model 
Ford fc Rasioll2006ft . iFaber et al.1 (|2005l ) finds that these 
violent passages might not destroy the planets, even if 
mass loss occurs. 



Previous studies have investigated gas free sys- 
tems wit h two planets around a cent ral sta r ex- 
tensively (iRasio fc Fordl [19961 : iFord et al.l l200lL [200l 
iFord fc Ra sio 2007) and have also begun to investigate 
the dynamics of two planets in the presence of a gas disk 
by using simplified prescri ptions for dissipative effects 
(|Moorhead fc Adams! |2005[ ). The observed eccentricity 
distribut ion is not easily reproduced by two equal-mass 
planets (jFord et al.1 l2001h . However, strong scattering 
of two unequal-mass planets could explain t he observed 
eccentricities of most ob served exoplanets (|Ford et al.l 
[200l IFord fc Rasidl200l . 

A system with three planets is qualitatively different 
than one with two planets. In two-planet systems, there 
is a sharp boundary between initial conditions that are 
provably Hill stable and i nitial conditions that quickly 
lead to a close encounter ( Gladmanl 1 1 993f ) . Moreover, 
for two Jupiter-mass planets in close to circular orbits, 
this boundary lies where the ratio of orbital periods is 
less than 2:1. If these planets formed further apart, then 
a slow and smooth migration could lead to systems be- 
comin g trapped in a 2:1 M MR bef ore triggering an insta - 
bility dLee fc Pealell200l but see iSandor fc Kle^l2f)M) . 
These stability properties are in sharp contrast to those 
of systems with three or more planets, for which there 
is no sharp stability boundary. These systems can be- 
come unstable even for much wider initial spacings and 
the t imescale to the first c lose encounter can be very 
long ([Chambers et al.1 H996) . Even if all pairs of adja- 
cent planets in the system are stable according to the 
two-planet criterion, the combined system can evolve in 
a chaotic (but apparently bounded) manner for an arbi- 
trarily long time period before instability sets in. This 
timescale to instability can easily exceed other timescales 
of interest here such as those for the formation of giant 
planets, orbital migration, or dispersal of the gas disk. 
Therefore, the stability properties of three-planet sys- 
tems can be studied with long-term orbital integrations 
without the need to implement any additional physics. 

A pioneering study of the stability and final or- 
bit al properties of three-plane t syste ms was performed 
by iMarzari" fc Weidenschillind (|2002t hereafter MW02). 
They explored the basic nature of instabilities arising in 
systems with three giant planets around a central solar- 
like star, and they determined the final orbital properties 
after one planet is ejected. This study used highly ide- 
alized initial conditions with three equal-mass planets or 
with one arbitrary mass distribution (middle planet twice 
as massive as the other two) . It was also computationally 
limited to a rather small ensemble of systems. 

Our two main goals in this new work were to extend 
the work of MW02 by using more realistic initial con- 
ditions (see fj2j) and to perform a much more extensive 
numerical survey in order to fully characterize the sta- 
tistical properties of outcomes for unstable three-planet 
systems. Given the chaotic nature of the dynamics, one 
needs to integrate many independent systems to char- 
acterize the statistical properties of the final planetary 
orbits (see a detailed analysis for the dependence of var- 
ious statistics on the sample size in Appendix (X[ see 
also lAdams fc Lau ghlin 2003). Each system has to be 
integrated for a long time, so that it reaches the orbit- 
crossing unstable phase and later evolves into a new, sta- 
ble configuration. Given the rapid increase in computing 
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power, we are now able to perform significantly more 
and longer integrations to obtain much better statistical 
results than was possible just a few years ago. 

In addition, we also present the results of new simu- 
lations for systems of three giant planets still embedded 
in a residual gas disk. Here our goal is different: we 
focus on the transition from gas-dominated to gas-free 
systems, in the hope of better justifying the (gas-free) 
initial conditions adopted in the first part of our work. 
However, implementing the physics of planet-disk inter- 
actions implies a considerably higher computational cost, 
preventing us from doing a complete statistical study of 
outcomes at this point. 

2. GAS-FREE SYSTEMS WITH THREE GIANT PLANETS 

In this section we consider systems with three unequal- 
mass giant planets orbiting a central star of mass 1 M@ at 
distances of several AU. The planets interact with each 
other through gravity and physical collisions only. We 
first present our assumptions and initial conditions, with 
particular emphasis on realistic mass distributions for the 
planets, and then we describe our numerical results and 
their implications. 

2.1. Initial Orbits 

For all systems the initial semimajor axis of the closest 
planet is always set at a i =3 AU. The other two planets 
are placed using the spacing law introduced by MW02, 

a-i+i = a t + KR H ,i, (1) 

where Rh,% is the Hill radius of the i th planet and we 
set K = 4.4 for all runs in this section. These choices 
are somewhat arbitrary, but are guided by the follow- 
ing considerations. F or a solar-mass centr al star, the ice 
line is around 3AU (|Kokubo fc Ida 2002) and it is dif- 
ficult to form giant planets closer to the star (see e.g., 
iKokubo fc Mai 12002ft. Although inward type II migra- 
tion fsee e.g.. lGoldreich fc Tremaineill980D can bring the 
giant planets closer to the star, we avoid putting the 
planets initially very close to the star since very small 
initial semi-major axis will lead to predominantly colli- 
sional outcomes. Furthermore, we would like to minimize 
the computing time, which leads us to consider closely 
spaced systems, while avoiding MMRs. 

Since our simulations in this section do not include the 
effects of gas, they are not intended to model the early 
phases of planet formation. Instead, at t = 0, we begin 
integrating fully formed planetary systems with a disk 
sufficiently depleted that the planets are free to interact 
with each other without significant dissipation from the 
disk. In §2.41 we show that the time until instability 
within a particular set of initial conditions does not af- 
fect the statistical properties of final outcomes. Indeed, 
we expect that the chaotic dynamics, both before and 
after the first close encounter, results in the distribu- 
tion of final outcomes being independent of the instabil- 
ity timescale. This justifies our choice of a very com- 
pact initial configuration with short instability timescale 
(~ 10 4 yr), which minimizes the computational cost. See 
MW02 and Appendix [B] for further discussion of the de- 
pendence of the instability timescale on K . 

Initial eccentricities are drawn from a uniform distribu- 
tion between — 0.1, and orbital inclinations are drawn 
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from a uniform distribution between 0° — 10° (with re- 
spect to the initial orbital plane of the innermost planet). 
To make sure that we could discern any inclination- 
dependent effects we used a slightly broader range of 
inclinations than seen in our solar system. However, in 
§2.61 we show numerically that the choice of initial incli- 
nations does not affect the distribution of final inclina- 
tions significantly. All initial phase angles are assigned 
random values between 0°-360°. 

2.2. Planetary Mass Distributions 

Our current understanding of planet formation remains 
full of uncertainties and no single prescription can claim 
to predict a correct planet mass distribution. For this 
reason, we consider three different prescriptions to con- 
struct plausible initial mass distributions for Jupiter-like 
planets. In all cases we adopt the standard core-accretion 
paradigm and we closely f ollow the simpl e plan et for- 
mation model described in Ko kubo fc Id a (2002 s here- 
after KI02). Planet masses depend on the distance of 
the planet from the central star through the gas surface 
density profile of the protoplanetary disk. 

2.2.1. Mass Distribution 1 

In this prescription, we first assign the planetary core 
masses M coro assuming a uniform distribution between 
1-10 M®. We assume that the cores accrete all gas within 
4 Hill radii (KI02) to reach a total mass M at a semima- 
jor axis a given by 

M = 27raAE gas + M core , (2) 

where A = is the feeding zone of the planet core and 
m is the Hill radius of the planet core, given by 

/ 1 M \ 1/3 

rff = U^rJ a (3) 

Here M* is the mass of the central star and a is the 
orbital radius of the core (assumed to be on a circular 
orbit). The gas surface density in the disk is given by 

_ 3 

S gas = / g Si (t^u) 2 S cm ~ 2 > ( 4 ) 

where the coefficient f g = 240 is the assumed gas-to- 
dust ratio (taken from KI02), Ei is the surface density 
at 1 AU, and the exponent comes from the minimum- 
mass disk model. We use Si = 10 in this case, which 
is a little higher than the minimum-mass Solar nebula 
value of 7. The choice of Ei is somewhat arbitrary and 
motivated to produce roughly Jupiter-mass planets. The 
initial masses of the planets obtained with this procedure 
are between about 0.4 Mj and 1.2 Mj. 

For this mass distribution we performed a set of 1000 
independent dynamical runs. 

2.2.2. Mass Distribution 2 

This is a slight refinement on the previous case in which 
we adopt an alternative prescription for accretion that 
takes explicitly into account the growing mass of the ini- 
tial planetary core. The initial core masses are chosen 
as in t|2.2.1l but the final mass of each planet is now de- 
termined using the following equations. Assuming that 
an infinitesimal mass dm accreted by a planet of mass m 
decreases the disk density by dE, we can write 

dm = —2nanHrHdY:, (5) 
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where tih is the number of Hill radii over which the mass 
is accreted. The final mass M of a planet at a distance 
a from the star, starting with a core mass M COTe can be 
obtained by integrating Eq. [5] as follows, 



~ 1/3 dm = -2Trn H a 2 



1 



3m: 



1/3 



<2E, 



(6) 



Si 



where Ej is the initial disk surface mass density. Solving 
Eq. [6] and replacing Ej from Eq. [4] we find 



M 



M* 1/3 3 4 /3 



f— 



-3/2 



2/3 



3/2 



(7) 



Here we use n# = 8, i.e., we assume that the core ac- 
cretes all mass within 4 Hill radii on either side. We use 
the same values for f g and Si from £12.2.11 

For this mass distribution we have integrated a smaller 
set of 224 systems. 

2.2.3. Mass Distribution 3 

We expect that the final distributions of different or- 
bital properties may vary significantly with different ini- 
tial mass distribution. To further test this mass depen- 
dence, we created a third set of systems with a broader 
planet mass distribution. Here we assign planetary 
masses exactly as in §2.2.11 but the initial core masses are 
chosen differently. We sample M core from a distribution 

1/5 

of masses between 1-100 M ffi uniform in M C ore, while as- 
suming again that these cores accrete all gas within 8 
Hill radii. The exponent in the core mass distribution 
and the surface density at 1 AU, Ei = 15, are chosen 
somewhat arbitrarily with the goal to obtain an initial 
mass distribution that peaks around a Jupiter mass but 
with a tail extending up to several Jupiter masses. The 
choices above produce initial masses spanning about an 
order of magnitude, in the range 0.4 Mj~AMj. More- 
over, the distribution for higher-mass planets resembles 
the mass distribution of observed exoplanets (see §2.8[) . 

For this mass distribution we have integrated a set of 
500 systems. 

2.3. Numerical Integrations 

We integrate each system for 10 7 yr, which is 2 x 10 6 
times the initial period of the closest planet (Ti^), and 
typically much longer than the timescale for the onset of 
instability. We u se the hybrid integrator of MERCURY6 . 2 
(Chambers 1999) and integrate the orbits symplectically 
while there is no close encounter, with a time-step of 10 
days, but switching to a Bulirsch-Stoer (BS) integration 
as soon as two planets have a close approach (defined to 
be closer than 3 Hill radii). Runs with poor energy con- 
servation (\AE/E\ > 0.001) with the hybrid integrator 
are repeated using the BS integrator throughout with the 
same \AE/E\ tolerance. This happens in ~ 30% of all 
runs, but our conclusions are not affected even if we re- 
ject these systems. We find that, in all systems, at least 
one planet is eventually ejected. Note that, for three- 
planet systems, following an ejection the remaining two 
planets may or may not be dynamically unstable. There- 
fore, we do not stop the integration following an ejection. 
Instead, we continue all integrations for two planets un- 
til a fixed stopping time of 10 7 yr. For systems with two 
remaining planets we check for Hill stability using the 
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Fig. 2. — Time evolution of semi-major axes and eccentricities for 
two randomly chosen typical simulations. The solid (black), dotted 
(red), and dashed (green) lines show the orbital elements for the 
initially closest (ai, ei), middle (122, e2), and furthest (0,3, e-i) 
planets. The top pair of panels show a realization where the first 
planet is ejected at ~ 4.1 X 10 4 T\ j, and the integration concludes 
with two planets in provably stable orbits. The semi-major axes 
for both P2 and P3 remain constant and the eccentricities oscillate 
stably on a secular timescale. The bottom pair shows another 
realization where P3 collides with Pi at ~ 4.2 X 10 3 T\ i; e2 keeps 
increasing until, a little before 10 5 T\,i Pi gets ejected, leaving a 
single planet in the system. Since a single orbit is always stable 
we stop the integration following this ejection. Numbers in the 
subscript represent the positional sequence of the planets starting 
from the star and letters "i" and "f" mean initial and final values, 
respectively, in all plots. 



known semi-analytic criterion (|Gladmanlll993f ). In our 
simulations about 9% of systems were not provably Hill 
stable at the integration stopping time. We discard those 
from our analysis. When a single planet remains (follow- 
ing a second ejection or when a collision took place), the 
integration is of course stopped immediately. 

For systems with two remaining dynamically stable 
planets the orbits can still evolve on a secular timescale 
(typically ~ 10 5 -10 6 yr for our simulated systems) 
much larger than the dyn amical i nstability timescale 
(I Adams k, Lau ghlin 20063; see also iMurrav fc DermottJ 
120001 Ch 71. We study these systems with two remain- 
ing stable planets by integrating the secular perturba- 
tion equations for a fu rther 10 9 yr with the analytical 
formalism developed in lFord et alj j2000). Note that the 
more standard formulatio n for Solar system dynamics 
(|Murrav fc DermottJ |2000T ) is not appropriate for these 
planetary systems because a significant fraction present 
very high eccentricities and inclinations. We find that 
for most of our simulated systems our chosen integra- 
tion stopping time effectively sampled the full parameter 
space (see detailed discussion in N2.9p . 

We treat collisions between planets in the following 
simple way ( "sticky-sphere" approximation) . A collision 
is assumed to happen when the distance between two 
planets becomes less than the sum of their physical radii. 
We assume Jupiter's density (1.33 gem -3 ) for all plan- 
ets when determining the radius from the mass. After 
a collision the two planets are replaced by a single one 
conserving mass and linear momentum. Because we ac- 
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Fig. 3. — Cumulative distributions showing initial and final ec- 
centricities of the planets. Top and bottom panels show the initial 
and final cumulative eccentricity distributions, respectively. In the 
top panel solid (black), dotted (red) and dashed (blue) lines rep- 
resent the closest, middle, and furthest planets, respectively. They 
are on top of each other because the initial eccentricity distribution 
is the same for all of the planets. In the bottom panel solid (black) 
and dotted (red) lines represent the final inner and outer planets, 
respectively. The dashed (blue) line shows all remaining planets in 
final stable orbits. 

count for collisions, our results are not strictly scale free. 
However, we find that collisions are relatively rare for our 
choice of initial conditions, so we still present all results 
with lengths scaled to and times scaled to T\^. 

Since Mass Distribution 1 corresponds to our largest 
set of runs, we first show our results from this set in 
detail in the following subsections ( ^2.41 - t j2.7[) . Results 
for the other two sets are summarized in H2.8I 

2.4. Overview of Results 

In Fig. [2] we show a couple of randomly selected, rep- 
resentative examples of the dynamical evolution of these 
systems, showing both chaotic phases and stable final 
configurations. Note the order-of-magnitude difference 
in timescale to first orbit crossing, illustrating the broad 
range of instability timescales (see also the discussion 
in Appendix [B]). We find that strong scattering be- 
tween planets increases the eccentricities very efficiently 
(Fig. [3]). The median of the eccentricity distribution for 
the final inner planets is 0.4. The median eccentricity 
for the final outer planets is 0.37, that for all simulated 
planets in their final stable orbits, is ~ 0.38. 

We compare our results with the observed eccentricity 
distribution of detected exoplanets in Fig. 2J For a more 
meaningful comparison we restrict our attention to ob- 
served planets with masses greater than 0.4 Mj, similar 
to the lower mass cut-off in our simulated systems. We 
also place an upper limit on the semimajor axis at 10 AU 
for the simulated final planet population to address the 
observational selection effects against discovering plan- 
ets with large orbital periods. Similarly, since planets 
close to the central star can be affected by additional 
physics beyond the scope of this study (e.g., tides, g en- 
eral relativistic effects; see I Adams fc L aughlin 2006b) we 
also omit observed close-in planets with semimajor axes 




Fig. 4. — Comparison between the simulated and observed ex- 
oplanet populations. The solid (black) line shows the cumulative 
distribution of the eccentricities of the remaining planets in their 
final stable orbits. The dashed (red) line is that for the observed 
population. For this comparison we employ a lower mass cut-off of 
0.4 Mj on the observed population addressing the fact that we do 
not have lower mass planets in our simulations. We also consider 
only the simulated planets that are finally within 10 AU from the 
star to address the fact that in the observed population we do not 
have planets further out. We also employ a lower semi-major axis 
cut-off of 0.1 ai t i on the observed population. 
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Fig. 5. — Cumulative frequency plots of semi-major axes of the 
initial (top panel) and final (bottom panel) planets. Vertical solid 
(black), dotted (red) and dashed (blue) lines show the initial val- 
ues in the top panel. These are vertical lines because the initial 
semi-major axes of the closest, middle and outer planets do not 
have a spread. Solid (black) and dotted (red) curves in the bot- 
tom panel show the final inner and outer planets' semi-major axes, 
respectively. 

below 0.1 ax,i- 

As seen in Fig. [3] and Fig. HJ our simulations slightly 
overestimate the eccentricities of the planetary orbits. 
However, the slopes of the cumulative eccentricity dis- 
tributions at higher eccentricity values are similar. In a 
realistic planetary system, there might be damping ef- 



6 Chatterjee, 




Fig. 6. — Final semi-major axis versus eccentricity plot. All 
lengths are scaled by the initial closest planet semi-major axis (here 
ai^i = 3.0AU). Black solid circles and red open stars represent 
the final inner and outer planets, respectively. Solid lines show 
different constant periapse lines with values O.f and 0.5. Note the 
high eccentricities and the close approaches towards the central 
star. The empty wedge shaped region in the a-e plane at high 
eccentricities is due to the requirement for orbital stability. 



TABLE 1 
Comparison of Eccentricity 
Distributions 





D 


P 


Mass distribution 1 


0.113 


0.15 


Mass distribution 2 


0.171 


0.01 


Mass distribution 3 


0.087 


0.32 



a For each mass distribution, we 
compare the final eccentricity dis- 
tribution of the simulated popula- 
tion with the observ ed exoplanet 
population (Figs. [4] 1181 ■ Using 
the kstwo function in Numerical 
Recipes, wc calculate the two- 
sample Kolmogorv-Smirnov statis- 
tic, D, and the corresponding prob- 
ability, P. In each case, the high 
value of P indicates that we can not 
reject the null hypothesis that both 
samples were drawn from the same 
population. 

fects from lingering gas, dust or planetesimals in a pro- 
toplanetary disk. While our simplified models already 
come close to matching the eccentricity distribution of 
observed planets, including damping may further im- 
prove this agreement. To be more quantitative, we per- 
form a Kolmogorov-Smirnov (KS) test and find that we 
cannot rule out the null hypothesis (that the two pop- 
ulations are drawn from the same distribution) at the 
85% level (Tabled]). In Sgm we will show that a broader 
initial distribution of planet masses results in an even 
better match to the observed eccentricity distribution. 

The top and bottom panels in Fig. [5] show the cumu- 
lative distributions of the initial vs final semi-major axes 
for the planets. The planet that is closest to the star ini- 
tially may not remain closest at the end of the dynam- 
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Fig. 7. — Left panel: cumulative frequency plots for the final 
eccentricities of the two subgroups, Group 1 (solid line) and Group 
2 (dashed line). Right panel: cumulative frequency plots for the 
final semi-major axes of the two subgroups, Group 1 (solid line) 
and Group 2 (dashed line). KS statistics results, D and P are also 
quoted for each of the above. D and P are as defined in Table [T] 

ical evolution. In fact, all three planets, independent 
of their initial positions, have roughly equal probabil- 
ity of becoming the innermost planet in the final stable 
configuration when the planet masses are not very dif- 
ferent. In 20% of the final stable systems, we find a 
single planet around the central star, two planets having 
been lost from the system either through some combi- 
nation of collisions and dynamical ejection. The other 
systems have two giant planets remaining in stable or- 
bits. We find that the planets in the outer orbits show a 
tendency for higher eccentricities correlating with larger 
semi-major axes (Fig.[6j) . We now know that many of 
the current obser ved exoplanets may have other planets 
in distant orbits (|Wright et alj 12007). From our results 
we expect that planets scattered into very distant bound 
orbits will have higher eccentricities. Long-term radial 
velocity monitoring should be able to test this predic- 
tion. 

Next, we investigate to what extent the final orbital 
properties depend on the instability timescale (cquiv- 
alently, on how closely packed the initial configuration 
was). For each system we integrated in §2.2. \\ we noted 
the first time when the semi-major axis of any one of 
the planets in the system changed by at least 10%. We 
use this as a measure of the dynamical instability growth 
timescale. Then, we divide the set into two subgroups, 
based on whether this growth time was below (Group 1) 
or above (Group 2) its median value (so 50% of the in- 
tegrated systems are in each group) . Fig. [7] compares 
the final eccentricity and semi-major axis distributions 
between the two groups. We find that the distributions 
are indistinguishable, demonstrating that the final (ob- 
servable) orbital properties are not sensitive to when ex- 
actly a particular system became dynamically unstable, 
as long as the dynamics was sufficiently active (ensur- 
ing that close encouters occur) and avoiding initial con- 
ditions so closely packed that physical collisions would 
become dominant. This result is hardly surprising since 
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Fig. 8. — Cumulative histogram of the pericenter distance of 
the initial (top panel) and final (bottom panel) planets bound to 
the star. In the top panel the solid, dotted and dashed lines show 
the pericenter distributions of the initial closest, middle and the 
furthest planets, respectively. In the bottom panel the solid and 
dotted lines show the same for the final stable inner and the outer 
planets with their semi-major axes less than 10 AU. The dashed 
magenta line shows the pericenter distribution of the observed ex- 
oplanet population for comparison purposes. 

we expect the chaotic evolution to efficiently erase any 
memory of the initial orbital parameters. Our results 
can therefore be taken as representative of the dynam- 
ical outcome for analogous systems with an even larger 
initial spacing between planets (but avoiding mean mo- 
tion resonances; see Appendix |B|) . In practice, perform- 
ing a large number of numerical integrations for these 
more widely spaced initial configurations would be pro- 
hibitively expensive (see Appendix [B| . 

2.5. Hot Jupiters from Planet-Planet Scattering 

We find that a significant fraction of systems emerge 
with planets in orbits having very small periastron dis- 
tances. Fig. [6] shows the final positions of the planets 
that are still bound to the central star in the a — e plane. 
The solid lines represent different constant pericenter dis- 
tances. Note that the planets show weak correlations 
between the eccentricity and the semi-major axis. For 
the inner planets, a lower semi-major axis tends to im- 
ply higher eccentricity, while the outer planets show an 
opposite trend. The final inner and outer planets form 
two clearly separated clusters of points in the a — e plane 
due to stability considerations. 

Fig. [5] shows the cumulative distribution of the pe- 
riastron distances of the final bound planets around 
the star. For the sake of comparison, we also show 
the pericenter distribution of all observed exoplanets 
in Fig. [8] We see that 10% of the systems harbor 
planets with periapse distances < 0.05(11^, whereas, 
a few (~ 2%) harbor planets with periapse distances 
< O.Olai^. Since we do not include tidal effects, we can- 
not compare this quantitatively with the observed pop- 
ulation. However, this is consistent with the ~ 5% of 
observed planets with semi-major axes within 0.03 AU. 
If the initial semi-major axes are sufficiently small tidal 
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Fig. 9. — Cumulative distribution showing initial and final or- 
bital inclinations of the planets with respect to the initial invari- 
able plane. In the top panel the solid (black) and the dotted (red) 
lines represent the initial and final RMS inclination distributions of 
the planet orbits with respect to the initial invariable plane. In the 
middle panel solid (black), dotted (red) and dashed (blue) lines rep- 
resent the closest, middle and furthest planets, respectively. The 
bottom panel shows the final orbital inclination distributions of the 
remaining planets in the system. The solid (black), and the dot- 
ted (red) lines represent the inner and outer planets, respectively. 
The dashed (blue) line represents the relative angles between the 
two remaining planetary orbits. Note that the final closer planets, 
which are the planets more easily observable in a planetary sys- 
tem, statistically have higher inclinations. Note, that the relative 
inclinations between the planetary orbits are also quite high. 

forces could then become important and a planet on 
a highly eccentric orbit could be circularized to pro- 
duce a hot Jup i ter dFord fc Rasio 2006; F aber et ap | 2005|: 
iRasio fc Fordl 119961; jW cidcnschillin g fe Marzaril Il99a 
iMarzari fc Weidenschillingl 12002) . However, recall that 
systems with much smaller values of a± t i would also lead 
to more physical collisions than in our simulations. More- 
over, a full numerical study of this scenario should in- 
clude tidal dissipation as part of the dynamical integra- 
tions, and possibly also include additiona l physics such 
as GR effects, etc. (jNagasawa et al.l[2008T l . 

2.6. Planets on High-inclination Orbits 

Since the star and planets get their angular momenta 
from the same source, planetary orbits are generally ex- 
pected to form in a coplanar disk perpendicular to the 
stellar spin axis. In Fig. [9l we compare the distributions 
of the final inclination angles. Here each angle reported 
is the absolute value of the orbital inclination measured 
with respect to the initial invariable plane, defined as 
the plane perpendicular to the initial total angular mo- 
mentum vector of the planetary orbits. Note that the 
direction of the initial total angular momentum can dif- 
fer from the direction of the total angular momentum of 
the bound planets at the end of a simulation, since plan- 
ets are frequently ejected from the system, carrying away 
angular momentum. 

Strong scattering between planets often increases incli- 
nations of the orbits, leading to a higher final RMS value 
of planet inclinations compared to the initial configura- 
tion (Fig. [9j top panel). In general, the inclinations tend 
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Fig. 10. — Initial RMS inclination vs final inclination of the inner- 
most planet. Note that the final closer planet orbital inclination is 
largely insensitive to the initial RMS inclination. 



to increase for all planets. The middle and bottom pan- 
els in Fig. [9] show the initial and final inclinations of the 
orbits of individual planets, respectively. The inclination 
of the final inner planet is typically larger than that of 
the final outer planet (Fig. [9) bottom panel). 

Our results show that strong planet-planet scattering 
can dramatically affect the coplanarity of some plane- 
tary systems (Fig. HI bottom panel). Since the timescale 
for tidal damping of incli nations is usually much greater 
than the age of the stars (|Winn et al.ll2005l ) , significantly 
increased inclinations could be found in some plane- 
tary systems that have gone through strong gravitational 
scattering phases in their lifetimes. Measuring a poor de- 
gree of alignment among the planetary orbits in multiple- 
planet systems, or between the angular momentum of one 
planet and the spin axis of its host star, could be used 
to identify systems that have undergone a particularly 
tumultuous dynamical history. 

If a system were initially assigned to a strictly copla- 
nar configuration, then angular momentum conservation 
dictates that it would remain coplanar always. How- 
ever, away from this trivial limit, we expect little corre- 
lation between the initial and final inclinations, given the 
chaotic nature of the dynamics. We test this hypothesis 
here by investigating the correlation between the initial 
and final inclinations of all planets in our simulations. 
We find that the final inclination of the inner planet 
indeed does not depend on the initial RMS inclination 
(Fig. [T0|) . We can quantify the amount of correlation 
between the initial RMS inclination and the final orbital 
inclination of the final inner planet using the bivariate 
correlation coefficient. The bivariate correlation coeffi- 
cient (r xy ) for two variables x and y, is given by the 
following equation, 



Cov(x,y) 
sd(x)sd(y) ' 



(8) 



where Cov(x, y) is the covariance of x and y, and sd(x) or 
sd(y) is the standard deviation of x or y. We find that the 
correlation coefficient between the initial RMS and the 
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Fig. 11. — Pericenter distance vs inclination of the final inner 
planets. The open dots show the final positions of the final inner 
planet in the pericenter-inclination plane. The filled disks (blue) 
and triangles (red) represent the mean orbital inclination of the in- 
ner planet and the final RMS inclinations, respectively. The means 
are obtained for bins of equal population (n(,j n = 50). We observe 
a weak anti-correlation between the pericenter and the inclination. 



final inner planet orbital inclinations is TmMS,iciose — 
0.05. The low value of r confirms that the high final 
inclinations are not merely a reflection of the initial con- 
ditions. As long as the planetary system is not strictly 
coplanar initially, strong planet-planet scattering can in- 
crease the orbital inclinations of some systems signifi- 
cantly. 

The final inclination of the inner planet, which is the 
most easily observable, shows a weak anti-correlation 
with the pericenter distance of its orbit (Fig. [TTj): lower 
pericenter orbits tend to have higher inclinations. The 
correlation coefficient in this case is r rp i C i ose = —0.13 
(Eq.EJ. 

For our solar system, the angle between the spin axis 
of the Sun and the invariable plane is ~ 6°. The an- 
gle between the stellar rotation axis and the orbital 
angular momentum of a transiting planet (A) can be 
constrained via the Rossiter-McLaughlin effect. Obser- 
vations have measured A sin i for fi ve systems (Win n 
2006b): —4.4° ±1.4° for HD 209458b (I Winn et al.l l2005h. 
-1.4°±1.1° for HP 189733b (IWinn et al.l l2006fl. 11°±15° 
for HD 149026b (IWqlf et al.l l2007ft. 30°±21° for TrES-lb 
(|Narita et apl2007aD, and, mos t recently, 62° ± 25° for 
HP 17156b (jNarita et alll2007bh . Our study implies that 
planetary systems with a tumultuous dynamical history 
will sometimes show a large A. Therefore, we look for- 
ward to precise measurements of A for many planetary 
systems to determine the fraction of planets among the 
exoplanet population with a significant inclination. In 
particular HP 17156b is very interesting in this regard, 
since the potentially high A together with the high ec- 
centricity (e = 0.67) strongly indicates a dynamical scat- 
tering origin for this planet. Measurements of A would 
be particularly interesting for the massive short-period 
planets (m > Mj), the very-short period giant plan- 
ets (P < 2.5 d), or the eccentric short period planets, 
since these planets might have a different formation his- 
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tory than the more common short-period planets with 
m ~ 0.5M,/ in nearly circular orbits. 

2.7. Mean Motion Resonances 

The radial-velocity planet population currently in- 
cludes 20 multi-planet systems and at least 5 of those sys- 
tems are in MMR (4 appear to be in a 2:1 MMR). MMRs 
can have strong effects on the dynamical evolution and 
stability of planetary systems. The 2:1 MMR is particu- 
larly interesting given the proximity of the two orbits and 
the increased possibility for close encounters that could 
result in strong gravitational scattering between t he two 
planets (]Sandor fc Klevll2lM Isihidor et al.ll200l . 

It is widely believed that MMRs between two or more 
planets in a planetary system arise naturally from mi- 
gration. Convergent migration in a dissipative disk can 
lead to resonant ca pture into a stabl e MMR, particu- 
larly the 2:1 MMR (jLee fc Peaiell2002t ). Simulations in- 
cluding an empirical dissipative force show that plan- 
etary orbits predominantly get trapped in 2:1 MMR 
(jMoorhead fc Adamdl2005T : iNagasawa et alJl2008f ). 

While we regard differential migration as a natural way 
to trap planets into MMRs, we did explore the possibility 
of trapping two planets into 2:1 MMR using only the mu- 
tual gravitational perturbations and without any damp- 
ing. We certainly expect this to be more difficult than 
with damping. Finding even a few systems trapped in 
MMR without any dissipation would be both surprising 
and interesting. In a three-planet system it is possible 
that one planet acts as a source or sink of energy to let 
the other two planets dynamically evolve into or out of 
a MMR. If pure dynamical trapping into MMRs were ef- 
ficient, then this would open up interesting possibilities. 
For one, it does not require a common disk origin, as is 
a requirement for the migratory origin of MMRs. Ad- 
ditionally, this mechanism could operate in a planetary 
system at a much later time after the protoplanetary disk 
has been dissipated. 

To look for possible 2:1 MMR candidates, we isolate 
systems that have two remaining planets with their final 
periods close to a 2:1 ratio. Then we calculate the two 
resonance angles 9\ and 82 over the full time of their 
dynamical evolution. Here the two resonance angles are 
given by 

01,2 =^l-2^2+roi, 2 , (9) 

where fa and fa are the mean longitudes of the inner and 
outer planets and w± and W2 are the longitudes of perias- 
tron for the inner and outer planets, respectively. When 
the planets are not in a MMR, 9%$ circulate through 2ir. 
When trapped in a MMR, the angles librate around two 
values (Lee 2004) . Finally, we check whether the periodic 
ratio and libration of the resonant angles are long lived 
or just a transient stage in their dynamical evolution. 

We find one system where two planets are clearly 
caught into a 2:1 MMR (Fig. [TJ]). The top two pan- 
els show the time evolution of the resonant arguments 9\ 
and 9i- The two resonant angles go from the circulating 
phase to the librating phase at around 1.88 x lO 6 ! 7 ^. 
The two bottom panels show the evolution of the semi- 
major axes and the eccentricities of the two planets in 
MMR. Note that the semi-major axes are nearly con- 
stant and the eccentricities oscillate stably. Since there 
is no damping in the system, the somewhat large libra- 
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Fig. 12. — Time evolution plots for the two resonance angles $1 
and 02, the semi-major axes and the eccentricities of the planets. 
From top to bottom the panels show the time evolutions of 0\, 62, 
semi-major axes and eccentricities, respectively. The time axis is 
in units of the initial orbital period of the initially closest planet 
(Tn). For the panels showing semi-major axes and eccentricity, 
the solid (red) and dotted (blue) lines show the evolutions of the 
two planets that enter a 2:1 MMR. Note that a little before 1.88 X 
10 6 Ti ti both 0i and 62 start librating. 

tion amplitude of the resonant angles is to be expected. 
In principle, the presence of even a little damping (due to 
some residual gas or dust in the disk) might reduce the 
amplitudes of libration and eccentricity oscillations for 
systems such as this one. A case like the one illustrated 
in Fig. [12] is clearly not a typical outcome of purely dy- 
namical evolution. We found a few other systems (~ 1%) 
showing similar librations of 9\^ at different times dur- 
ing their dynamical evolution, but only for a brief phase 
never exceeding ~ 1Q A T\^. However, if our simulations 
had included even some weak dissipation, the frequency 
of such resonances might have increased significantly. We 
encourage future investigation of this possibility. 

2.8. Mass Dependences 

Our simulations show the effects of mass segregation, 
as heavier planets preferentially end with smaller semi- 
major axes. This trend can be easily seen by compar- 
ing the initial and final mass distributions of the planets 
in Fig. [131 The mass distribution clearly shifts towards 
higher mass values in the final inner planet mass his- 
togram, whereas, the outer planet mass more closely re- 
flects the initial mass distribution (compare the top and 
bottom panels of Fig. IT3|) . We do not find a strong effect 
of mass on eccentricity but we note that collisions tend to 
reduce the fraction of highly eccentric systems (Fig. [14]) . 
The collision products can be seen in the cluster around 
and above 1.5 Mj. We find no other significant mass de- 
pendent effect in the final orbital parameters for our set 
of runs using Mass Distribution 1 . 

We now describe briefly the results obtained with the 
two alternative initial mass distributions for the three 
planets. Fig. [15] shows correlation between semi-major 
axis and mass for both Mass Distribution 1 ( §2.2.ip and 
Mass Distribution 2 ( ij2.2.2p . Somewhat surprisingly, for 
Mass Distribution 2, we find no significant differences 
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Fig. 13. — Initial and final mass distributions of the closest, mid- 
dle and furthest planets. The top (bottom) panel shows the initial 
(final) mass distributions. Solid (black), dotted (red), and dashed 
(blue) lines in the top panel represent the initial mass histograms 
of the closest, middle, and furthest planets. Solid (black) and dot- 
ted (red) lines in the bottom panel represent the mass histograms 
of the final inner and outer planets, respectively. One planet is 
ejected in each of our simulations. Note that the histogram for 
the inner planet masses shifts towards higher values in the bottom 
panel, which indicates that the higher mass planets preferentially 
become the inner planet in the final stable configuration of the 
planetary systems. 




Fig. 14. — Mass vs eccentricity of the final stable planets. The 
circles (black) and the triangles (blue) represent the final inner 
and the outer planets. Planets with masses l.GMj are collision 
products. The collision planets tend to have lower eccentricities. 
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Fig. 15. — Mass vs semi-major axis of the final remaining stable 
planets. The disks (black) and the open stars (red) represent the 
final inner and the outer planets for Mass Distribution 1. The open 
squares (blue) and the open circles (green) represent the same, 
respectively, for Mass Distribution 2. Note, higher mass planets 
remain close to their initial positions. 




Mass (Mj 

Fig. 16. — Same as Fig. 1131 but with a different and wider initial 
mass distribution than the previous one (Mass Distribution 3). 
The initial mass distribution has a high number of Jovian mass 
planets as in the previous mass distribution, however, in this case 
the distribution has a tail towards higher masses. The higher end 
in the initial mass spectra in this case mimics the minimum mass 
(msini) spectrum of the observed exoplanets. Note that the mass 
segregation effect is more prominent here than in Fig. 1131 The 
dot-dash (green) line shows the m sin i distribution of the observed 
exoplanets in both panels for comparison. 



from the results obtained with the much simpler prescrip- 
tion of Mass Distribution 1 . This is possibly because in 
both Mass Distributions 1 and 2, the mass range and dis- 
tribution are similar (Mass Distribution 2 is only shifted 
towards slightly higher values). 

For this reason we also studied a third choice of mass 
distribution, Mass Distribution 3 f tj2.2.3[) . with a much 
larger range of planetary masses, enabling us to observe 



mass-dependent effects more clearly. For example, we 
now see that the tendency for higher mass planets pref- 
erentially to become the final inner planets (Fig. [T6|) is 
more prominent than in our other simulations. Similarly, 
the effect of a mass distribution on the final eccentric- 
ities of the remaining planets is more prominent with 
this broader mass distribution. The higher-mass planets 
preferentially excite the eccentricities of the lower-mass 
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Fig. 17. — Same as Fig. [3] but using Mass Distribution 3. Note Fig. 18. — Same as Fig. [4] but using Mass Distribution 3. The 

that overall the eccentricities are lowered using the broader mass simulated eccentricities match much better with the observed in 

distribution. this case compared to those using Mass Distribution 1. 



counterparts, often to the point of ejection. This effec- 
tively reduces the overall eccentricities of the final stable 
orbits (Fig. I17p . The median value of the final inner orbit 
eccentricities is 0.24, and that for the outer orbit is 0.23 
in this case. The final cumulative distribution of eccen- 
tricities matches the observations even more closely with 
Mass Distribution 3 than with Mass Distribution 1 or 2 
(Fig.[l8l Table[T|). We employ similar selection criteria as 
described in H2.5I The final semi-major axis distribution 
is statistically indistinguishable from the one obtained 
with Mass Distribution 1. We also clearly see that the 
lower-mass planets get scattered around preferentially 
while the heavier counterparts do not move much and 
stay mostly near their initial positions (Fig. flT))) . This is 
in accord with the observation that close-in planets are 
often of low er mass than planets with moderate s emi- 
major axes (|Cumming et al.ll2lM iNaef et alJ l2005). At 
present, the correlation between planet mass and orbital 
period for radial- velocity planets is consistent with a pop- 
ulation of systems where the less massive planets have 
been scattered inwards. We predict that planet searches 
sensitive to longer-period planets will eventually find a 
population of sub- Jupiter planets that have been scat- 
tered outwards. Furthermore, our simulations predict 
a negative correlation between mass and orbital period 
among such long-period planets, if they are launched into 
their current orbits via strong gravitational scattering. 
We find that mass and eccentricity have a weak anti- 
correlation (Fig. [20]) . We do not find any systems with 
two planets trapped in 2:1 MMR for this case. 

2.9. Secular evolution 

It is known from numerous previous studies that sec- 
ular perturbations of one planet on another in a multi- 
planet system can modify the planets' orbital properties 
on a timescale much longer than the relevant dynami- 
cal (orbital, or strong dyn amical in stability) timescales 
(I Adams fc Laughlinl [2006a; see also iMurrav &: Dermottl 
2000). Since secular timescales can be orders of magni- 
tude longer than the orbital timescales, one might obtain 
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Fig. 19. — Same as Fig. 1151 but using Mass Distribution 3. Mass 
dependent effects on final semi-major axes of the planets is much 
more prominent in this case. 



results biased towards the initial part of the oscillations 
if at least a full secular period is not sampled properly. 
Fig. shows a dramatic example where the eccentrici- 
ties of both planets and the relative inclinations between 
the planetary orbits oscillate secularly with a very long 
period (~ 100 Myr) compared to the orbital timescale 
and the observed eccentricities and inclinations can be 
very different from what would be expected right after 
dynamical stabilization of the system. Hence, any study 
of orbital properties of planets after dynamical interac- 
tions should also worry about the secular evolution of 
the orbital properties that follows the orders of magni- 
tude quicker dynamical phase. Nevertheless, we should 
point out that in our simulated systems this is not typ- 
ical. For most cases the secular time period is typically 
~ 10 5 - 10 6 yr. For our simulated systems containing 
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Fig. 20. — Same as Fig. 1141 but using the broader distribution 
of initial planet masses (mass distribution 3). There seems to be 
a weak anti-correlation between the mass and the eccentricities of 
the planets. 




Fig. 22. — Scatter plot of eccentricity after secular evolution for 



10 yr vs eccentricity after the integration stopping time ( i]2.3 
Left and right panels show the inner and outer planets, respec- 
tively. Note that the eccentricities for the inner orbits change more 
significantly than for the outer ones. 
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Fig. 21. — Secular time evolution of the relative inclination (top 
panel) and the eccentricities (bottom panel) of a system with two 
dynamically stable planets, t = for this is the end of dynamical 
integration f i|2.3l l. In the bottom panel the curve with a smaller 
oscillation amplitude (red) shows the evolution of the outer planet 
eccentricity and the other (black) shows that of the inner planet. 
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Fig. 23. — Left panel: Cumulative distributions of eccentricities 
of th e inner planets before and after secular evolution described in 
£|2.9I Solid (black) and dashed (blue) lines show the distributions 
before and after secular evolution, respectively. Right panel: Same 
as left panel, for the outer planet. KS test results for both pairs of 
distributions are also shown in the plot. 



two provably stable planets at the end of our common 
integration stopping time (10 7 yr) we study the evolu- 
tion of the eccentricities for a further 10 9 yr to confirm 
that the orbital properties at the end of our integration 
correctly represent the true final distribution. 

To evaluate the secular evolution of these plan- 
ets, w^e_jase_theoctupole-order formalism presented 
by iFord et all (|200Ct ). Note that the more stan- 
dard formulation in term s of the Laplace coefficients 
(|Murrav fc Dermottl |2000| ) is not appropriate for these 
planetary systems because a significant fraction of these 
systems contain orbits with very high eccentricities and 
inclinations. 



We find that indeed individual eccentricities of these 
planetary orbits can change significantly. Fig. [22] shows 
a scatter plot of the final eccentricities after secular evo- 
lution for 10 9 yr as a function of the eccentricities after 
our integration stopping time for both planets. It is clear 
that the individual eccentricities can change significantly, 
especially, for the inner planet. However, the overall dis- 
tribution does not change significantly from the distribu- 
tion obtained right after our integration stopping time in 
§2.31 Fig. l23l shows that the eccentricity distributions be- 
fore and after secular evolution for the outer planet, in 
particular, are statistically identical. For the inner plan- 
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Fig. 24. — Left panel: Cumulative distributions of relative in- 
clinations between the plan etary orbits before and after secular 
evolution described in £|2.9I Solid (black) and dashed (blue) lines 
show the distributions before and after secular evolution, respec- 
tively. KS test results for the two distributions (before and after 
secular evolution) are also shown in the plot. 



ets we find that, after secular evolution, there is a little 
overabundance of very high eccentricity (e 0.8) orbits 
(Fig. r25|) . In order to quantify the likeness of the two dis- 
tributions before and after secular evolution, we perform 
KS tests for both the inner and outer planet eccentric- 
ity distribution. We find that we cannot rule out the 
null hypothesis (that the distributions before and after 
secular evolution are drawn from the same distribution) 
at 62% and 1% significance level for the inner and outer 
planetary orbits, respectively. The very low values of the 
significance level along with the large ensemble essen- 
tially means that the two distributions are very similar. 
We perform the same test with the relative inclination 
of the planetary orbits in the subset of our systems with 
two dynamically stable remaining planets (Fig. l24l) . For 
these distributions the significance level for KS test with 
the same null hypothesis is 27%. This confirms that our 
choice of integration stopping time already sampled the 
full parameter space for the secular evolution. 

3. EFFECTS OF A RESIDUAL GAS DISK 

In the previous section we considered the dynami- 
cal evolution of three-planet systems with fully formed 
planets on initially near-circular orbits and no gas disk. 
Implicit assumptions are that sufficiently massive disks 
damp planetary eccentricities, and that residual gas disks 
dissipate quickly enough to allow the later chaotic evo- 
lution of planetary systems. Here, we will verify these 
assumptions by simulating three-planet systems within 
residual gas disks. 

3.1. Photoevaporation 

The final stage of disk dissipation remains poorly 
understood. Since viscous evolution alone cannot ex- 
plain the observed rapid disp ersal of disks (~ 10 5 yr; 
see e.g., iSimon fc Pratoiri995f) . some other mechanism 
must be responsible for removing a res idual disk . The 
most likely is photoevaporation (e.g., IShu et ail 119931 : 



Fig. 25. — Surface mass densities for runs with DISK1-9 
(black line). Also shown is the standard surface mass density 
£ = 10 3 (r/AC/)- 3 / 2 (yellow line). The surface mass densities for 
DISK 1—9 are obtained by evolving a disk with the standard surface 
mass density under disk's viscosity parameter of a = 5 X 10 — 3 for 
3, 5, 10, 13, 15, 17, 20, 25, and 30 Myr. 



iHollenbach et al.l I1994D . iClarke et all (|2001l ) proposed 
that, once the viscous accretion rate drops to a level 
comparable to the wind mass loss rate, photoevapora- 
tion takes over the disk evolution. When this limit is 
reached, surface layers of the disk beyond the gravita- 
tional radius (R g — GM/c 2 s ), where the sound speed 
c s exceeds the disk's escape speed, starts removing disk 
mass faster than it is being replenished by viscous evo- 
lution. As a result, the disk is divided into inner and 
outer parts: the inner disk drains onto the central star 
on a short viscous timescale, while t he outer disk evap- 
orates on longer tim es cales (e.g., IClarke et al.l 120011 : 
I Alexander et al.ll2006h . I Alexander et alj (|2006h showed 
that the disk clearing by this mechanism takes about 
10 5 yr, which is comparable to the observed dissipation 
time. 

The viscous evolution time at semi-major axis a is de- 
fined as 

t v is{a) = — -7- — - (10) 



A^disk (a) 
Miisk - 3ttz/E , 



(11) 



where v and £ are the viscosity and surface mass density, 
respectively. 

On the other hand, the photoevaporation time at a is 



iphoto(fl) 



M disk (< a) 



(12) 



' wind (^0 

w here the wind mass loss rate for an optically thick disk 
is (|Clarke et al.ll2001l ) 

lO^J I Mo J Meyr_ ' 

(13) 

and for an optically thin disk ((Alexander et a l. 2006) 



M wind = 4.4x10" 



M wind = 9.68 x 10" 
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TABLE 2 
Disk models 





Disk Mass 


JJisk Age (1U yr ) 


tdamp (10 u yr) 


DISK1 


3.7 Mj 


0.3 


0.2 for 17/24 


DISK2 


2.2 M j 


0.5 


0.4 for 14/23 


DISK3 


0.22 Mj 


1.0 


4 for 17/28 


DISK4 


23.5 M 9 


1.3 


7 for 5/27 


DISK5 


11.8 M s 


1.5 




DISK6 


4.7 M© 


1.7 




DISK7 


1.4 M© 


2.0 




DISKS 


0.19 M© 


2.5 




DISK9 


0.02 M s 


.3.0 





a The disk masses of the 9 different disk models used in g3] 
The disk age is the time until a typical disk with a = 5 X lCr 
will reach the corresponding total mass. The eccentricity 
damping time is the time to reduce the planetary eccentric- 
ities from above 0.1 to below 0.1, and obtained for systems 
which did not go through mergers. 



«in \ 
3 AU/ 



1/2 



\ a out 



M e yT-\U) 



Here $ is the ionizing flux from the central star, h is 
the pressure scale height of the disk, ai n and a out are the 
inner and outer disk radii. 

Photoevaporation becomes effective when i V j S > iphoto- 
For typical disks, this corresponds to a disk mass of a 
few Jupiter masses. When a disk mass drops below this 
critical value, planets are likely to become dynamically 
unstable if the photoevaporation time is shorter than the 
dynamical instability growth time (iphoto < idyn)' In this 
section we will investigate this further by simulating 3- 
planet systems with various disk masses. 

3.2. Numerical Method and Assumptions 

For this study we use a hybrid A-body and 1-D gas dy- 
namics code to follow the evolution of three-planet sys- 
tems for several different disk masses. Our hybrid code in 
its current form combines an existing A-body integrator 
with a 1- D implementatio n of a viscous, nearly Keplerian 
gas dis k (lThommesll2005D. The A-body code is based on 
SyMBA (|Duncan et al.lll99cf ). It is fast for near-Keplerian 
systems, requiring only ~ 10 timesteps per shortest orbit, 
while undergoing no secular growth in energy error. In 
addition, it makes use of an adaptive timestep to resolve 
close encounters between pairs of bodies. 

The gas disk is divided into radial bins, each of 
which represents an annulus whose properties (sur- 
face density, viscosity, temperature, etc.) are az- 
imuthally and ve rtically averaged, follow ing the gen- 
eral approach of iLin fe Papaloizoul (|1986f ). Arbitrary 
viscosities can be specified through a st andard a- 
parametrization (|Shakura fc Svunvaevl [T9731. Though 
this disk is explicitly 1-D, the vertical and azimuthal 
structures are implicitly included in the model. For 
the former, a scale height is assigned to every annu- 
lus. The latter is key to the planet-disk interactions, 
which result from the raising of azimuthally asymmet- 
ric structure (spiral density waves) in the disk by the 
planet. This effect is added in the form of the tor que 
density presc ription of iGoldreich fc Tremaind (|1980f ). as 
modified by IWardl (|1997| ). which describes the disk- 
planet angular momentum exchange taking place as 
waves are launched. Planetary eccentricities are damped 
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FlG. 26. — Instability growth timescale for 30 planetary sys- 
tems with DISK1-9 (from right to left). The range of dynami- 
cal instability time is relatively independent of disk masses. Data 
points with an arrow indicate the number of systems which did 
not go through dynamical instability within our simulation time 
(~ 10 7 yr). Diagonal lines are disk clearing times by photoevapo- 
ration for optically thin disks (dark lines) and optically thick disks 
(light lines). Solid and dashed lines are the disk clearing time mea- 
sured at the gravitational radius R g and 10 AU respectively. The 
horizontal line shows the viscous evolution timescale of the disk 
at R g . Photoevaporation is expected to take over disk evolution 
once t v i scoua (R g ) ~ iphoto(^s)- Planetary systems are likely to go 
through a chaotic evolution as shown in gas free systems ([[2} when 
tphoto < *dyn (i- e - above t photo lines in the figure). 



on timescales as in IWardl (H99l and I Artvmowic3 (fl993T) . 
Since we do not take account of the saturation of coro- 
tation resonances, which could lead to the eccentric- 
ity excitation by Lindblad reson ances (jGoldreich b, SarH 
l2003t iMoorhead fc Adams 2008) , the eccentricity damp- 
ing considered here is an upper limit. 

3.3. Results: Onset of Dynamical Instability 

For initial conditions, we randomly choose 30 three- 
planet systems from the set using Mass Distribution 1 
in ij2.2.11 and study their orbital evolution within 9 dif- 
ferent disk masses. The surface mass density profiles of 
these disks are shown in Fig. These are obtained by 
evolving a minimum mass solar nebula disk model with a 
viscosity parameter a — 0.005 for various times (without 
planets). Disk properties are summarized in Tabled and 
we will refer to our models as DISK 1-9 from here on. We 
assume that each of these disks is inviscid for dynamical 
runs with planets, meaning that type II planet migration 
is not taken into account. However, this should not af- 
fect our results significantly since even the most massive 
disk (DISK1) contains only 3.7Mj, which is comparable 
to the planetary masses used in our simulations. Most of 
our disks are therefore too small to affect planet migra- 
tion. 

The dynamical instability is commonly characterized 
by the orbital crossings of planets. Fig. [2^1 shows the first 
orbital crossing time of each system for each disk mass. 
Diagonal lines are disk clearing timescales by photoevap- 
oration for optically thick disks (orange/grey lines), and 
optically thin disks (blue/black lines). Also plotted is the 
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viscous evolution times at the gravitational radius. This 
figure indicates that photoevaporation takes over disk 
evolution for disks with a few to several Jupiter masses, 
depending on the photoevaporation models. 

It appears that the range of first orbital crossing time 
idyn is relatively independent of disk masses, and around 
~ 10 — 10 4 yr. Note, however, that the number of systems 
going through orbital crossings decreases for larger disk 
masses. Excluding mergers, 19/24, 14/28, and 4/28 sys- 
tems for DISK1, 2, and 3, respectively do not experience 
any orbital crossings during the simulation time (10 7 yr), 
while all systems with lighter disk masses go through at 
least one crossing within the run time. This is concor- 
dant to expectations that planets become dynamically 
unstable more readily in a less massive gas disk and the 
same planetary system that remained stable in a suffi- 
ciently massive disk can become unstable once the disk 
dissipates. 

Apart from td yn , the eccentricity damping timescale 
(idamp) is very important to know, since tdamp determines 
whether the disk can damp the eccentricities back to near 
zero after one (or more) orbit-crossing episode(s), before 
the disk is depleted. We define tdamp as the time taken to 
damp the planetary eccentricities from e > 0.1 to e < 0.1. 
For 30 different systems for 4 disk masses (DISK1-4) we 
find the median tdamp to be 2 x 10 5 , 4 x 10 5 , 4 x 10 6 , and 
7 x 10 6 yr, respectively. For less massive disks we do not 
find significant damping. While the evolution of disks is 
dominated by viscous evolution and tdamp < t v is (DISK1 
and 2), planetary orbits are expected to remain nearly 
circular since after an instability there is enough time to 
damp the eccentricities before the disk is depleted. The 
nature of evolution can change drastically once photoe- 
vaporation dominates the disk evolution and starts de- 
pleting the disk more rapidly. We find that most systems 
reach at least one orbit-crossing episode for the least mas- 
sive disks (DISK8 and 9) since td yn > tphoto- Some plan- 
etary systems in more massive disks have tdyn < tphoto 
(Fig. |21)|) . For these more massive disks eccentricities 
excited via planet-planet interaction may be damped if 
tdamp < ^photo- However, since the median tdamp is longer 
than tphoto for these disks, planetary eccentricities ex- 
cited via planet-planet interaction do not have time to 
be damped before the gas disk is depleted, once photoe- 
vaporation is efficient. 

In summary, we expect that planetary systems will re- 
main stable with nearly circular orbits while the planets 
are embedded in a sufficiently massive disk. Even if there 
is an occasional orbital crossing or merger, the eccentric- 
ities and inclinations will rapidly damp in such a disk, so 
that the system returns to nearly circular orbits. Then 
eccentricities will evolve more freely once photoevapora- 
tion takes over the disk evolution, and the disk clearing 
time becomes short compared to the instability growth 
time (tdyn > tphoto)- Even when planets become unsta- 
ble before the disk is completely depleted (tdyn < tphoto), 
it is unlikely that their eccentricities are damped, since 
the eccentricity damping times of these disks tend to be 
longer than the disk dissipation time (tdamp > t p hoto)- 
Therefore, we expect that most planetary systems be- 
come dynamically unstable when a gas disk dissipates. 
This further justifies our initial conditions in Sj2j In a 
future paper we will further investigate the evolution 



of multiple-planet syste ms within an evolving gas disk 
(jMatsumura et al.ll2008h . 

4. COMPARISON WITH PREVIOUS STUDIES 

The previous work most similar to ours was the pio- 
neering study by MW02 on (gas-free) three-planet sys- 
tems. MW02 also studied the orbital properties of plan- 
etary systems following a dynamically active phase of 
their evolution. However, their study was computation- 
ally limited and their systems were rather idealized in 
terms of assumed planetary masses and initial orbits. 
Our results are in good qualitative agreement with those 
of MW02. For example, they showed for the first time 
with three-planet systems how scattering can produce 
large eccentricities. However, our more realistic and gen- 
eralized initial conditions enable us to explore a larger 
parameter space and to study in more detail the most 
interesting phenomena such as the generation of large, 
potentially observable inclinations. We also find that the 
final stable planets can be scattered at even smaller semi- 
major axes than they predicted. Since these very low 
semi-major axes planets are in the tail of the distribu- 
tion, it is expected that a simulation of a smaller sample 
size will miss some of them (see Appendix |A|. Moreover, 
our much larger simulated sets and improved statistics 
on dynamical outcomes allow us to better compare our 
theoretical predictions to observations (see Appendix [A|. 

In addition to the orbital properties of remaining plan- 
ets, MW02 also presented a stability timescale analy- 
sis for planetary systems with three giant planets. In 
verifying these results, we realized the importance of 
this study, especially for our choice of initial spacing, 
and we therefore decided to perform a much more de- 
tailed timescale analysis, with significant improvements 
over MW02 made possible by the dramatically increased 
speed of present-day computers. The results of this anal- 
y sis are presented in A p pendi x iBl 

iMoorhead fc Adamsl (|2005l ) studied in detail two- 
planet systems with an empirical dissipation arising from 
a residual disk. They found that, even in initially well 
separated two-planet systems, migration can bring the 
planets close enough for dynamical instability. In their 
study they accounted for a disk outside both planets with 
their empirical formula, whereas, we immerse the three 
planets in a protoplanetary disk with varying disk masses 
C - Another major difference between their study and 
ours is the number of planets considered. The dynami- 
cal evolution of two-planet systems can be very different 
from that of systems with three or more planets (see 
jLj. Keeping these differences in mind, we compare key 
points between the two studies. For example, for suffi- 
ciently massive disks we find that the eccentricity damp- 
ing timescale is less than the disk dissipation timescale. 
However, as the disk mass is diminished, the timescale 
for eccentricity damping and the number of unstable sys- 
tems increases. They also find that scattering fills up the 
a-e plane for the inner planet orbit. Due to the setup 
of their initial conditions and the dynamical limitations 
of two-planet system they do not find planets with large 
orbital periods, normally produced by strong scattering 
between planets. They also stop integrating after 1 Myr 
or when the system has only one planet left. One should 
remember that in cases where two planets are remaining, 
the planetary properties can still change either through 
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dynamical scattering (see discussion in £|2.3[) or even for 
dynamically stable systems, through long term secular 
perturbations ( for a detailed discu s sion s ee N2.9|) . 

More recently I Juric fc Tremainel (|2007f ) perform an in- 
teresting s t udy a s an extension of the pioneering work by 
iLin fc Idal ()1997l ). They study the dynamical evolution 
of generic iV-planet systems with N > 3 and a wide range 
of initial conditions. Although their three-planet systems 
were dynamically inactive as a result of their choice of 
initial separations and integration stopping time, their 
other runs with higher TV bear very relevant results for 
our study. In particular, they find a similar final ec- 
centricity distribution, suggesting that this distribution 
may be universal. One of the most interesting results 
in their study is that the final number of surviving plan- 
ets following a dynamically active phase is almost always 
2-3, independent of the initial number. Since these sys- 
tems are chaotic, the properties of any planetary system 
emerging out of a dynamically active phase will have lit- 
tle memory of the initial number of planets or the exact 
initial conditions (also see discussion in i j2.4l and fc|2.6|) . 
One can imagine a situation where a system started with 
N > 3 and, followed by many collisions and ejections, 
reaches a stage with N = 3. If dissipation circularizes 
orbits after each ejection or collision, then such a system 
could reach a state similiar to the initial conditions for 
our three-planet simulations. Thus, our results may be 
r epresentative of e v en mo re generic multi-planet systems. 

iNagasawa et al.l (|2008l ) study the dynamics of three 
equal-mass planets including dynamical tides. Since they 
can apply tides while the three-planet dynamical scatter- 
ing phase is still active, they find increased efficiency to 
tidally isolate planets that would otherwise still actively 
take part in three-planet scattering. In this study we did 
not include tides. However, we find that 8% of these sys- 
tems have one planet accreted onto the star. We also find 
that ~ 40% of these systems contained at least one planet 
which, during the three-planet violent scattering phase, 
reached a pericenter distance within 0.01 AU of the star. 
In most of the systems these close-in planets end up be- 
ing ejected, while, some of them collide with the star or 
another planet. None of these systems remain stable at 
the end of the run. The addition of tides during these 
scattering phases may stabilize some systems by isolating 
the close-in planet from the other planets dynamically. 
Of course the circularization process needs to be very 
efficient and quick so that the planet gets circularized 
and decoupled from the others before it can be e j ected . 
We should also point out that Na gasawa et alj (2008) 
study equal-mass planets. We find that the dynamical 
evolution of planetary systems with unequal masses is 
very different than for equal mass systems (see discus- 
sion in £12.21 fc E|2 .8[) . In particular, the lower-mass plan- 
ets preferentially get scattered inwards or outwards while 
the heavier counterparts remain near their initial posi- 
tions throughout the whole evolution (e.g., see Fig. [T5|). 
One should also remember that the tidal circularization 
timescale depends on the mass and ra dius of these plan- 
ets (e.g., Ilvanov fc Papalol zou 2004). This can affect 
the efficiency of tidal circularization for high eccentric- 
ity planets in their setup. At present their study actually 
produces too many (30%) "hot" planets compared to the 
current observed population of ~ 5% within 0.03 AU (at 
0.03 AU the tidal circularization timescale is ~ 10 6 yr for 



Jovian planets, see Nagasawa et al. 2008). Note that the 
selection biases of radial velocity surveys can only reduce 
the fraction of hot Jovian planets in the future. It will be 
interesting to see results of similar studies with a more 
realistic mass distribution. 

5. SUMMARY AND CONCLUSIONS 

We have studied in detail how the orbital properties 
evolve through strong gravitational scattering between 
multiple giant planets in a planetary system containing 
three giant planets around a solar mass star. We per- 
form a detailed study for gas free generic planetary sys- 
tems. We focus on the final orbital properties of the 
planets that remain bound to the central star in sta- 
ble orbits after chaotic evolution due to strong mutual 
interactions, followed by a prolonged secular evolution 
(~ 10 9 yr) when two planets remain after the scattering 
phase. We perform the experiments with realistic plan- 
etary systems containing 3 giant planets (fj2j). In all of 
these systems at least one planet is eventually ejected 
before reaching a stable configuration. This supports 
models of planet formation that predict planetary sys- 
tems initially form several closely spaced planets, but 
instabilities reduce the number of planets until the sta- 
bility timescale exceeds the age of the planetary system. 
In 20% of the cases, two planets are lost through ejec- 
tions or collisions leaving the system with only one giant 
planet. Thus, the planet scattering model predicts the 
existence of many systems with a single eccentric giant 
planet, as well as many free floating planets (depend- 
ing on how many planets are formed before the planet 
scattering phase of evolution) . 

We find that strong gravitational scattering between 
giant planets can naturally create high-eccentricity or- 
bits. The exact distribution of eccentricities for the final 
remaining planets in stable orbits depends on the choice 
and range of the initial mass distribution. When the ini- 
tial mass distribution spans a broad range of masses, the 
less massive planets typically start to acquire larger ec- 
centricities. However, these planets with highly excited 
orbits arc often ejected, reducing the overall eccentricities 
of the remaining dynamically stable planetary orbits. Al- 
though the first two sets of our models (with a narrower 
range of initial planet masses) predict eccentric planets to 
be slightly more common than observed (Fig. 2]), a wider 
initial mass distribution can result in remarkable simi- 
larity with the observed distribution (Fig. [T8|) . Recently, 
a s imilar trend in ecc e ntricit ies was found independently 
by I Juric fc TremaTn c (2007) for generic dynamically ac- 
tive multi-planet systems independent of the details of 
the initial conditions or the initial number of planets. 

We conclude that planet-planet scattering could eas- 
ily account for the observed distribution of eccentricities 
exceeding 0.2. However, our simulations slightly under- 
produce systems with eccentricities less than 0.2. This 
may suggest that some observed systems are affected by 
late stage giant collisions. Alternatively, the presence of 
a residual gas or planetesimal disk could lead to eccen- 
tricity damping. We find this latter explanation partic- 
ularly attractive given the ob served correlation between 
planet mass and eccentricity (|Butler et al.ll2006D . While 
our simulations suggest that high eccentricities are most 
common among less massive giant planets, the known 
population of extrasolar planets suggest that high ec- 
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centricit ies are more comm on among the more massive 
planets (|Ford fe Ra sio 2007). This apparent discrepancy 
could be resolved if a modest disk often remains after the 
final major planet-planet scattering event. Less massive 
planets would be more strongly affected by the remain- 
ing disk, so their eccentricities could be damped, while 
more massive planets would typically be immune to ec- 
centricity damping. 

We find that it is possible to scatter some planets into 
orbits with low perihelion distances (Fig. [H]). Approxi- 
mately 10% of the systems obtain perihelion distances 
less than 0.05 ai^, whereas, a fewer fraction (~ 2%) can 
reach within O.Olai.i. If the initial semi-major axes are 
small enough, then strong gravitational scattering could 
result in planet orbits with sufficiently small perihelion 
distances, such that tidal effects could circularize their 
orbits at small orbital distances. 

We find that the inclination distribution of such plan- 
ets could be significantly broadened. If we assume that 
the angular momentum of the host star is aligned with 
that of the initial orbital angular momentum of the plan- 
ets, then measurements of A (the angle between stel- 
lar spin axis and planet's orbital angular momentum) 
should typically be small in absence of perturbations 
from other planetary or stellar companions f £|2.6[) . We 
find that strong gravitational scattering between the gi- 
ant planets can naturally increase the inclinations of 
the final planetary orbits with respect to the initial to- 
tal orbital angular momentum plane (Fig. [9|). Since 
the timescale to tidally align the stellar spin and the 
planetary angular momentu m is much great e r than the 
age of the star (~ 10 12 yr; lGreenberdfl97l lHutlll98(l 
I Winn et "ail 120051 ). inclinations excited by planet-planet 
scattering after the disk had dispersed could be main- 
tained for the entire stellar lifetime. Observations of 
a hot- Jupiter with a significantly non-zero A would be 
suggestive of previous planet-planet scattering. How- 
ever, caution woul d be necessary if the star had a binary 
stellar companion (IWu fc Murrav l l2003: iTakeda fc Rasiol 



I2005t iFabrvckv fc Tremaind 120071 ). On the other hand, 
observations of many hot-Jupiters with orbital angular 
momenta closely aligned with their stellar rotation axis 
would suggest a formation mechanism other than strong 
gravitational scattering followed by tidal circularization. 
Unfortunately, current observations measure this angle 
for only a few systems and some measurements have 
uncertainties comparable to the dispersion of inclina- 
tions found in our simulations. We encourage observers 
to improve both the number and precision of Rossiter- 
McLaughlin observations. 

We find that the relative inclinations between plane- 
tary orbits in the systems with two remaining planets 
in their final stable orbits also increase via planet-planet 
scattering. Future observations using astrometry or tran- 
sit timing could possibly measure relative inclinations 
between planetary orbits in multi-planet systems. Fur- 
thermore, we find that in ~ 20% of the systems hav- 
ing two giant planets in their final dynamically stable 
configurations, the relative inclination between the two 
planets is higher than 40°. For these systems it is possi- 
ble for the planets to g o through Kozai-type oscillations 
(Nagasawa ct al. 2008). Although effects of a debris disk 
on planetary dynamics and vice versa is beyond the scope 
of this study, the warped disk observed in (3 Pictoris could 



be one interesting example where inclined planetary or- 
bits and the debris d isk exchange torques, resulting in 
a war p ed debris disk (ISmith fc Terrildll984l; iHeap et all 
I2000D . iMouillet et ail (|1997f) suggests that the observed 
asymmetry in the debris disk can be explained by the 
presence of a planetary companion in an inclined orbit. 
Strong planetary scattering, as we find, can be a natu- 
ral way to create planetary orbits with large semi-major 
axes and highly inclined orbits. 

Less massive planets are more likely to be scattered 
far away from the site of their formation. Our simula- 
tions show that both the close-in or farther out planets 
should have lower mass than the planets with moderate 
semi-major axes (Figs. [131 [19")) . This trend can be veri- 
fied in future observations using adaptive optics to detect 
and image giant planets further out (4 0-100 AU) from 
the central star (|Lafreniere et al.ll2007f ). We find that a 
few percent of the simulated population have very high 
semi-major axes in the final stable configuration (e.g., 
Fig-EJ). Such giant planets are extremely unlikely to be 
created in situ, since the timescale fo r planet formation 
great ly exceeds the age of the star ()Veras fe Armitagel 
l2004f ) . Additionally, there is simply insufficient disk mass 
to form a giant planet at such large orbital distances 
dKokubo fc Idall200l Hda fc Linl I2004b llah . Strong scat- 
tering between planets in multi-planet systems can be 
a natural mechanism to create such long-period planets 
(a > 50 AU). Our simulations suggest that this popula- 
tion of high semi-major axis planets will have high eccen- 
tricities and inclinations (Fig. [6|) . Future planet searches 
using astrometry or direct detection can test these pre- 
dictions. 

We have also presented a preliminary study of the ef- 
fects of a residual gas disk on planetary dynamics ($3}. 
We compare the importance of dynamics for 9 different 
disk models with different disk surface densities, keep- 
ing the initial orbital properties of the embedded planets 
the same for all cases. We identify important timescales 
for the dynamical evolution of these systems. In par- 
ticular, we characterize the transitional stage of the dy- 
namical evolution from the stable, eccentricity damped 
phase, where the planets are embedded in a massive disk 
to the unstable free eccentricity evolution stage following 
disk depletion. We show that it is possible to understand 
the overall evolution after planets are fully formed as an 
interplay between four different timescales, namely: the 
viscous timescale (t V i s ) of the disk, the photo-evaporation 
timescale (t photo) of the disk, the dynamical instability 
timescale td yn for the planetary orbits, and the eccentric- 
ity damping timescale (tdamp) for the planetary orbits in 
a disk ( §3.3p . Our study clearly shows that planets will 
remain stable on nearly circular orbits while a sufficient 
amount of gas remains present in the disk, while with 
the same initial orbits without gas the system would be- 
come unstable. The boundary between these two differ- 
ent phases can be characterized by tdamp, t V is, and t p hoto- 
The unstable phase starts when gas mass is sufficiently 
depleted so that tdamp > t p hoto/vis- We find that the 
transition within a disk from a small to a large tdamp 
can be fairly quick once t p hoto < tvi a (Fig. [26)) . After 
photo-evaporation takes over the disk evolution, the sys- 
tem undergoes a quick transition. Until the critical mass 
for photo-evaporation is reached, planetary eccentricities 
remain close to zero independent of the disk mass and 
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previous dynamical history. Then, once the critical den- 
sity is reached, the system behaves as if it had started 
with near-circular initial planetary orbits in a gas-free 
environment. Thus, apart from highlighting the relative 
importance of these timescales for the evolution of plane- 
tary systems, our results justify typical initial conditions 
used in most studies of gas-free multi-planet systems, in- 
cluding our own (j2J). The initial properties for a gas- free 
system would be the orbital properties of the system as 



found at the boundary where tdamp > tphoto/visi m this 
context. 
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APPENDIX 

SAMPLE SIZE NEEDED TO ACCURATELY CHARACTERIZE THE ECCENTRICITY DISTRIBUTION 

Due to the chaotic nature of n-body integrations and finite precision of computer arithmetic, the long-term integration 
of an exact system is impossible. Instead, we must perform ensembles of n-body integrations and interpret the results 
in terms of the statistical properties of the outcomes of a set of similar n-body systems. Given that any study is based 
on a finite sample size, it is important to recognize the limitations on the precision of various statistics due to the 
limited number of integrations. In the context of this paper (and similar works) many simulations of similar planetary 
systems are performed to generate a random sample of outcomes, each of which can be compared to the observable 
properties, such as the eccentricity and the semi-major axis of extrasolar planetary systems. Here we present an 
analysis of precision of several statistics describing the eccentricity distribution as a function of the number of n-body 
integrations performed. Formally, our estimated precision is applicable only to our specific choice for the distribution 
of initial conditions. Since this and previous studies have shown that many variations of the planet-planet scattering 
model result in similar eccentricity distributions, we expect that our results can be applied to many similar studies to 
estimate the accuracy of various statistical properties. Our results should also give a quantitative way to decide the 
required number of simulations to estimate various statistical properties to a given precision, and we expect this will be 
of use to other researchers when formulating research plans. Since the eccentricities of the planetary systems are one 
of the most interesting properties, we focus our attention on statistics describing the distribution of final eccentricities. 
However, the basic idea can be applied to any statistic describing the masses or orbital properties of the simulated 
planetary systems. 

We will estimate the precision of several interesting statistics (the mean, the standard deviation, the 5th (P5) and 
the 95th (P95) percentiles of the distribution of final eccentricities). We estimate the "true" value for the underlying 
population based on estimates obtained making use of our full sample of N = 1515 simulations (see fj2]). We then 
estimate the same statistics based on m subsets of the full sample, where each subset is a random sample of n systems. 
Next, we compare the statistic estimated from each subset of simulations to the statistic based on the full population. 
We systematically change the sample size n for each subset, so as to explore how precisely we can estimate a given 
statistic as a function of the sample size. 

In Fig. [27j we show each estimate for a given statistic as a single tick mark. Each panel presents results for a 
different statistic: mean (left, top), standard deviation (left, bottom), P95 (right, top), and P5 (right, bottom) of the 
eccentricity distribution. For each of the above statistics, we show the mean plus and minus the standard deviation 
(blue short dashed curves) and the 5th and 95th percentiles (magenta long-dashed curves). Hence, the upper magenta 
curve of the upper left panel shows the P95 for the estimate of the mean eccentricity based on a sample of m — N/n 
estimates of the mean eccentricity each using a sample size of n. We find that for n = 50 the standard error in 
estimating the population mean from the sample can have a standard error of 37%, whereas, for n = 100 the precision 
improves to 5% (Fig. |27|) . In estimating the standard deviation of the eccentricity from the small subsamples, we 
find that with n — 100 the standard error in the estimation of the standard deviation of the eccentricity distribution 
is ~ 7% (Fig. |27|) . Although estimates for first two moments of the eccentricity distribution can come within 10% 
of the population moments based on only n ^ 100, a much larger sample size is required to characterize the tails of 
the eccentricity distribution accurately. For example, the errors in estimating the P5 and the P95 of the underlying 
eccentricity distribution using sample sizes of n — 100 are 30% and 7%, respectively. If the sample size is increased to 
n ~ 1000, then P5 and P95 can be estimated within 4% and 1% error, respectively, 

Both previous and future studies often generate predictions for the eccentricity distribution of planetary systems 
based on various theoretical models. For the sake of comparing the precision with which these studies estimated the 
predicted eccentricity distribution, we have used the above result to obtain an empirical relation between the standard 
deviation of the estimates of various summary statistics describing the eccentricity distribution and the number of 
simulations used to estimate the statistic (Fig. |2"8"]) , We expect that this relation will also be useful for planning 
future studies, where researchers will want to make a deliberate choice regarding the number of simulations and other 
simulation parameters such as length of integration time, number of particles, and inclusion of additional physics. 
We find that the standard deviation in the deviation of the estimated mean eccentricity from the population mean 
eccentricity decreases as a power law of sample size n with an index of —1.586 ± 0.004. The standard deviation 
estimating P5 and P95 decreases less steeply with n; here the power law indices are —0.58 ± 0.05 and 0.51 ± 0.03, 
respectively. (These empirical relations are valid only for n ^ 5.) For example, a study that uses 100 simulations would 
typically estimate the mean of the predicted eccentricity distribution to within ~ 0.008. However, a larger number 
of simulations becomes increasingly important for estimating the tails of the eccentricity distribution precisely. For 
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Fig. 27. — Top left panel: Each dash shows the mean eccentricities estimated from randomly chosen samples of sizes n, specified by the 
x-axis. The solid (red) line shows the mean eccentricity of the full sample. The short dashed (blue) lines show the standard deviation in 
the estimates of the mean eccentricity with respect to the full population mean for each sample size. The long dashed (magenta) lines show 
the 5th and the 95th percentiles of the estimates of the mean eccentricities. Bottom, left panel: The same for the estimates of the standard 
deviation for the eccentricities. Top right panel: The same for the estimates of the 95th percentiles for the distribution of eccentricities. 
Bottom right panel: The same for the estimates of the 5th percentiles for the distribution of eccentricities. 



example, an ensemble of 100 simulations typically estimates the P5 or the P95 with a precision of ~ 0.025 or ~ 0.045, 
respectively. 

STABILITY TIMESCALE 

According to the core accretion model of plane t formation, planets form in a protoplanetary disk separated by a 
small number of Hill radii away from each other (|Kokubo k. Idal [ l9 98, 2002). Hence, it is very interesting to have a 
good and statistically reliable investigation of the stability timescales as well as the distributions of the timescales as 
a function of the planet-planet distances in multiples (K) of their mutual Hill radii. A similar timescale study was 
also performed by I Chambers et al.l (|1996f) . However, their study covers a very different range of planetary masses. 
The large ensembles used by our study not only produce a better statistical characterization of these timescales as 
a function of K, they also enable us for the first time to show the actual nontrivial shapes of these distributions. 
The actual distributions of these timescales for a given K value are particularly interesting for anyone performing a 
similar study and trying to decide upon a reasonable initial planetary separation, since due to the broad range, the 
computational effort needed will be determined by the few unusually stable realizations rather than the more frequent 
ones where instability can grow orders of magnitude quicker. For better comparison with MW02 we put the planet 
closest to the star at 5AU and then determine the semi- major axes of the other two planets as follows, 

a%+\ = a% + KRff^i+i, (Bl) 

where K is the spacing measured in terms of Rn,i,i+i, the mutual Hill radius for the i th and i + 1 th planets, 

Mi + M i+ i \ 1/3 at + a i+ i 



'{-^mt 1 ) ^-t^' (B2) 

following their prescription. Here Mi is the mass of the i th planet, M* the mass of the central star, and a,i the semi- 
major axis of the i th planet. Note we use a different definition of Hill radius from that in §2.21 following MW02 for 
easier comparison. 

We integrate a number of 3-planet systems with different initial conditions: 1000 for K ^ 4.3, 500 for 4.3 < K ^ 5.0 
and 200 for K ^ 5.0. Apart from the large number of realizations, we use a more general distribution of the initial 
eccentricities and orbital inclinations as described in M2.21 

Fig. [29] shows the results as a function of K. The filled circles show the median and the vertical bars above and 
below represent ±34% around the median. Note that the vertical bars are not error bars, but they are representative 
of the actual distributions of the stability timescales. We also show the mean of each distribution to compare it with 
the median. In each case the mean overestimates the timescale and lies often outside the 34% bars around the median. 
Our results are consistent with the findings of MW02 qualitatively. We see similar trends near a MMR. However, we 
find that a simple linear fit as was tried by MW02 does not work well. A better empirical fit is given as follows. 

logiot m (K) = a + b x exp(cK), (B3) 
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Fig. 28. — Standard deviations in the estimates of the mean (top panel), the 5th percentile (bottom left panel), and the 95th percentile 
(bottom right panel) for the eccentricity distribution based on a randomly chosen sample of size n. The solid (red) lines show the empirical 
fit for the standard deviation of the estimation of the statistic of interest as a function of the sample size. The empirical formulae for 
the fitting lines are also shown at the top of each panel. Note that the fits are only valid for samples of size greater than 5. Moreover, 
the standard deviation in the estimates for very large n (and hence very small standard deviation) is limited by the precision of our main 
simulated data (and is the reason for the saturation seen in the standard deviation values near n = 1000). 

where a, b, and c are constants (henceforth, called Fit-1). The best fit values for a, b, and c (Table [3j can predict 
the median timescales with fractional error less than 10% away from MMR. We also tried to find a simpler linear fit 
t m ,imear(K) for our data away from MMR, following MW02, writing 

loijiotjn, linear 

{K) = a + bK, (B4) 

where a and b are fitting parameters (henceforth, called Fit-2). The best fit values for a and b are also given in Table [31 
We find that Fit-1 is much better than Fit-2. In particular, we find that linear fitting formula for instability timescale 
as a function of initial spacing (as suggested by MW02) is inaccurate by over three orders of magnitude for initial 
spacings such that the planets are beyond the 2:1 MMR. 

We also present the first study of the actual shapes of the timescale distributions. In particular, in cases of broad 
or skewed distributions, knowing only the median (or mean) timescale can not provide a complete description of the 
distribution of timescales to instability. We find that the shapes of the distributions of the timescales arc essentially 
the same for any K value away from a major MMR whereas near a major MMR the shape is very qualitatively different 
with a much slower decay above the median timescale (Figs. [30 l l3T j) . Both systems near and away from MMR show a 
similar exponential part in the stability timescale distribution. However, due to the MMR configuration, some of the 
systems enjoy increased stability manifested as a broader distribution to the higher time end. Note that the histograms 
of the timescale distributions are normalized such that J^. ruAlti = 1, where AZtj is the bin size in logarithm of time. 
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Fig. 29. — Instability growth timescale as a function of spacing parameter K. The filled circles are the medians of the distributions of 
stability timescales for the respective values of K . The bars show the ±34% range in the timescale. Note that they are not error bars, 
instead they represent the distribution of the stability timescales. The medians are shown in the plot because the distributions are skewed 
towards greater timescales. The open stars show the mean of the timescale distribution. Signatures of some major resonances can be 
noticed in the sudden dips in timescale. The arrow indicates that most systems with K = 5.5 are stable for at least 10 9 yr. The starting 
point of the arrow represents the shortest instability timescale among our simulations. The solid line (blue) and the da shed lin e (green) 
show the empirical best fit lines predicting the medians of the stability timescale distributions given the K values (Eqs. IB3IIB4I I. 



The normalized number distribution for times lower than the median timescale (henceforth denoted as ul) has an 
exponential shape; above the median timescale (henceforth denoted as ur) the number distribution has a linear decay 
for all K away from major MMRs. The fitting formulae for til and ur are given by 

n L = N L exp [(log 10 1 - log 10 t m {K))/t L ] , (B5) 

nR= N R - t R log w t. (B6) 

Here, Nl and Nr are the normalization constants for the peak amplitudes of the distributions, t m (K) is the median 
of the timescale distribution as a function of K , ti, and tR are fitting constants characterizing the exponential index 
and the slope of the two curves, respectively. The best-fit values for Nl, Nr, t^, and tR are listed in Table |4j For a 
given K value, the median timescale can be estimated using Eq. IB3l and then using the median timescale the shapes 
of the distributions can be obtained using Eqs. IB 5 1 and IB61 
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Fig. 30. — Histograms for the timescale distributions at two different K values both away from MMR. Note that times are shown in 
log scale. Each histogram shown here corresponds to 10 3 runs for that K value. The number distributions are normalized such that 
■ riiAlti = 1, where, Aitj is the bin-size in logarithm of time. This normalization essentially makes the area under each histogram 
normalized to 1. The solid histogram corresponds to K = 3.9 and the dashed histogram corresponds to K = 2.3. Both these K values are 
away from any major MMR. The two histograms have essentially the same shape. The dotted (blue and red) curves show the analytical 
fitting curves for timescale distributions at the left and the right sides of the mode of the distributions. For systems with stability timescales 
less than the medi an of th e distribution show an exponential shape, whereas, those with timescales higher than the median show a linear 
drop-off (see Eqs. [B5][B6]and Table B}. 



TABLE 3 







Fit: 


t vs K 






a 


b 


c 


Max Error(%) 


Fit 1 


1.07 


0.03 


1.10 


10 


Fit 2 


-1.74 


1.29 




50 



a The best fit values of the fitting parameters 
for the empirical fits for the median stability 
timescale of the systems as a function of their 
initial spacing parameter K. 
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Fig. 31. — Histograms for the timescale distributions near and away from an MMR. Each histogram corresponds to 10 3 runs for that K 
value. We follow the same normalization scheme as mentioned earlier. K = 3.3 is near the K value for a 3 : 2 commcnsurability between 
the periods of the first and the second as well as the second and the third planetary orbits (dashed line). K = 3.9 is away from MMR (solid 
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TABLE 4 
Fit: Timescale distribution 

( n L,R) 



K 


t L 


N L 


t R 


Nr 


2.0 


0.719 


0.714 


1.728 


3.716 


2.1 


0.669 


0.775 


1.138 


2.649 


2.2 


0.582 


0.886 


0.778 


1.971 


2.3 


0.595 


0.864 


0.836 


2.129 


2.4 


0.663 


0.867 


1.193 


2.816 


2.5 


0.530 


0.996 


1.559 


3.628 


2.6 


0.482 


1.153 


1.286 


3.065 


2.7 


0.589 


1.004 


1.112 


2.806 


2.8 


0.582 


1.002 


1.167 


3.01 


2.9 


0.713 


0.800 


0.584 


1.791 


3.0 


0.624 


1.715 


1.412 


3.141 


3.1 


0.810 


0.512 


0.095 


0.522 


3.2 


2.343 


0.206 


0.082 


0.484 


3.3 


1.447 


0.283 


0.203 


0.968 


3.4 


1.051 


0.462 


0.188 


0.889 


3.5 


0.533 


1.238 


0.180 


0.869 


3.6 


0.395 


3.700 


0.201 


0.929 


3.7 


0.811 


0.937 


0.143 


0.753 


3.8 


0.744 


1.085 


0.245 


1.140 


3.9 


1.211 


0.415 


0.522 


2.409 


4.0 


1.488 


0.353 


0.349 


1.811 


4.1 


1.399 


0.379 


0.160 


0.975 


4.3 


1.684 


0.346 


0.204 


1.292 


4.4 


2.650 


0.218 


0.127 


0.953 


4.8 


1.066 


4.910 


0.102 


0.710 


5.0 


0.767 


73.831 


0.209 


1.242 



a Best fit values for the fitting param- 
eters, Nl, t^, Nr, and tR predicting 
til, and ur for a given K and the me- 
dian of the stability timescale distri- 
bution t m (K). 
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